home *** CD-ROM | disk | FTP | other *** search
/ The Programmer Disk / The Programmer Disk (Microforum).iso / xpro / c2 / pro3 / tan.c < prev   
Text File  |  1987-06-20  |  6KB  |  282 lines
  1. /*                            tan.c
  2.  *
  3.  *    Circular tangent
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, tan();
  10.  *
  11.  * y = tan( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns the circular tangent of the radian argument x.
  18.  *
  19.  * Range reduction is modulo pi/4.  A rational function
  20.  *       x + x**3 P(x**2)/Q(x**2)
  21.  * is employed in the basic interval [0, pi/4].
  22.  *
  23.  *
  24.  *
  25.  * ACCURACY:
  26.  *
  27.  *                      Relative error:
  28.  * arithmetic   range      # trials      peak         rms
  29.  *    DEC      0,+1.07e9      44000     4.0e-17     1.0e-17
  30.  *    IEEE     0,10           30000     2.9e-16     8.0e-17
  31.  *    IEEE     0,2.147e9       6000     2.9e-16     8.3e-17
  32.  *
  33.  * ERROR MESSAGES:
  34.  *
  35.  *   message         condition          value returned
  36.  * tan total loss   x > 1.073741824e9     0.0
  37.  *
  38.  */
  39. /*                            cot.c
  40.  *
  41.  *    Circular cotangent
  42.  *
  43.  *
  44.  *
  45.  * SYNOPSIS:
  46.  *
  47.  * double x, y, cot();
  48.  *
  49.  * y = cot( x );
  50.  *
  51.  *
  52.  *
  53.  * DESCRIPTION:
  54.  *
  55.  * Returns the circular cotangent of the radian argument x.
  56.  *
  57.  * Range reduction is modulo pi/4.  A rational function
  58.  *       x + x**3 P(x**2)/Q(x**2)
  59.  * is employed in the basic interval [0, pi/4].
  60.  *
  61.  *
  62.  *
  63.  * ACCURACY:
  64.  *
  65.  * See tan().
  66.  *
  67.  *
  68.  * ERROR MESSAGES:
  69.  *
  70.  *   message         condition          value returned
  71.  * cot total loss   x > 1.073741824e9       0.0
  72.  * cot singularity  x = 0                  MAXNUM
  73.  *
  74.  */
  75.  
  76. /*
  77. Cephes Math Library Release 2.0:  April, 1987
  78. Copyright 1984, 1987 by Stephen L. Moshier
  79. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  80. */
  81.  
  82. #include "mconf.h"
  83.  
  84. #ifdef UNK
  85. static double P[] = {
  86. -1.30936939181383777646E4,
  87.  1.15351664838587416140E6,
  88. -1.79565251976484877988E7
  89. };
  90. static double Q[] = {
  91. /* 1.00000000000000000000E0,*/
  92.  1.36812963470692954678E4,
  93. -1.32089234440210967447E6,
  94.  2.50083801823357915839E7,
  95. -5.38695755929454629881E7
  96. };
  97. static double DP1 = 7.853981554508209228515625E-1;
  98. static double DP2 = 7.94662735614792836714E-9;
  99. static double DP3 = 3.06161699786838294307E-17;
  100. static double lossth = 1.073741824e9;
  101. #endif
  102.  
  103. #ifdef DEC
  104. static short P[] = {
  105. 0143514,0113306,0111171,0174674,
  106. 0045214,0147545,0027744,0167346,
  107. 0146210,0177526,0114514,0105660
  108. };
  109. static short Q[] = {
  110. /*0040200,0000000,0000000,0000000,*/
  111. 0043525,0142457,0072633,0025617,
  112. 0145241,0036742,0140525,0162256,
  113. 0046276,0146176,0013526,0143573,
  114. 0146515,0077401,0162762,0150607
  115. };
  116. /*  7.853981629014015197753906250000E-1 */
  117. static short P1[] = {0040111,0007732,0120000,0000000,};
  118. /*  4.960467869796758577649598009884E-10 */
  119. static short P2[] = {0030410,0055060,0100000,0000000,};
  120. /*  2.860594363054915898381331279295E-18 */
  121. static short P3[] = {0021523,0011431,0105056,0001560,};
  122. #define DP1 *(double *)P1
  123. #define DP2 *(double *)P2
  124. #define DP3 *(double *)P3
  125. static double lossth = 1.073741824e9;
  126. #endif
  127.  
  128. #ifdef IBMPC
  129. static short P[] = {
  130. 0x3f38,0xd24f,0x92d8,0xc0c9,
  131. 0x9ddd,0xa5fc,0x99ec,0x4131,
  132. 0x9176,0xd329,0x1fea,0xc171
  133. };
  134. static short Q[] = {
  135. /*0x0000,0x0000,0x0000,0x3ff0,*/
  136. 0x6572,0xeeb3,0xb8a5,0x40ca,
  137. 0xbc96,0x582a,0x27bc,0xc134,
  138. 0xd8ef,0xc2ea,0xd98f,0x4177,
  139. 0x5a31,0x3cbe,0xafe0,0xc189
  140. };
  141. /*
  142.   7.85398125648498535156E-1,
  143.   3.77489470793079817668E-8,
  144.   2.69515142907905952645E-15,
  145. */
  146. static short P1[] = {0x0000,0x4000,0x21fb,0x3fe9};
  147. static short P2[] = {0x0000,0x0000,0x442d,0x3e64};
  148. static short P3[] = {0x5170,0x98cc,0x4698,0x3ce8};
  149. #define DP1 *(double *)P1
  150. #define DP2 *(double *)P2
  151. #define DP3 *(double *)P3
  152. static double lossth = 1.073741824e9;
  153. #endif
  154.  
  155. #ifdef MIEEE
  156. static short P[] = {
  157. 0xc0c9,0x92d8,0xd24f,0x3f38,
  158. 0x4131,0x99ec,0xa5fc,0x9ddd,
  159. 0xc171,0x1fea,0xd329,0x9176
  160. };
  161. static short Q[] = {
  162. 0x40ca,0xb8a5,0xeeb3,0x6572,
  163. 0xc134,0x27bc,0x582a,0xbc96,
  164. 0x4177,0xd98f,0xc2ea,0xd8ef,
  165. 0xc189,0xafe0,0x3cbe,0x5a31
  166. };
  167. static short P1[] = {
  168. 0x3fe9,0x21fb,0x4000,0x0000
  169. };
  170. static short P2[] = {
  171. 0x3e64,0x442d,0x0000,0x0000
  172. };
  173. static short P3[] = {
  174. 0x3ce8,0x4698,0x98cc,0x5170,
  175. };
  176. #define DP1 *(double *)P1
  177. #define DP2 *(double *)P2
  178. #define DP3 *(double *)P3
  179. static double lossth = 1.073741824e9;
  180. #endif
  181.  
  182. extern double MAXNUM;
  183. extern double PIO4;
  184.  
  185. double tan(x)
  186. double x;
  187. {
  188. double tancot();
  189.  
  190. return( tancot(x,0) );
  191. }
  192.  
  193.  
  194. double cot(x)
  195. double x;
  196. {
  197. double tancot();
  198.  
  199. if( x == 0.0 )
  200.     {
  201.     mtherr( "cot", SING );
  202.     return( MAXNUM );
  203.     }
  204. return( tancot(x,1) );
  205. }
  206.  
  207.  
  208. static double tancot( xx, cotflg )
  209. double xx;
  210. int cotflg;
  211. {
  212. double x, y, z, zz;
  213. int j, sign;
  214. double polevl(), p1evl(), floor(), ldexp();
  215.  
  216. /* make argument positive but save the sign */
  217. if( xx < 0 )
  218.     {
  219.     x = -xx;
  220.     sign = -1;
  221.     }
  222. else
  223.     {
  224.     x = xx;
  225.     sign = 1;
  226.     }
  227.  
  228. if( x > lossth )
  229.     {
  230.     if( cotflg )
  231.         mtherr( "cot", TLOSS );
  232.     else
  233.         mtherr( "tan", TLOSS );
  234.     return(0.0);
  235.     }
  236.  
  237. /* compute x mod PIO4 */
  238. y = floor( x/PIO4 );
  239.  
  240. /* strip high bits of integer part */
  241. z = ldexp( y, -3 );
  242. z = floor(z);        /* integer part of y/8 */
  243. z = y - ldexp( z, 3 );  /* y - 16 * (y/16) */
  244.  
  245. /* integer and fractional part modulo one octant */
  246. j = z;
  247.  
  248. /* map zeros and singularities to origin */
  249. if( j & 1 )
  250.     {
  251.     j += 1;
  252.     y += 1.0;
  253.     }
  254.  
  255. z = ((x - y * DP1) - y * DP2) - y * DP3;
  256.  
  257. zz = z * z;
  258.  
  259. if( zz > 1.0e-14 )
  260.     y = z  +  z * (zz * polevl( zz, P, 2 )/p1evl(zz, Q, 4));
  261. else
  262.     y = z;
  263.     
  264. if( j & 2 )
  265.     {
  266.     if( cotflg )
  267.         y = -y;
  268.     else
  269.         y = -1.0/y;
  270.     }
  271. else
  272.     {
  273.     if( cotflg )
  274.         y = 1.0/y;
  275.     }
  276.  
  277. if( sign < 0 )
  278.     y = -y;
  279.  
  280. return( y );
  281. }
  282.